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Abstract We present a multiscale approach to simulate the impact of a solid object on 
a liquid surface: upon impact a thin liquid sheet is thrown upwards all around the rim 
of the impactor while in its wake a large surface cavity forms. Under the influence of 
hydrostatic pressure the cavity immediately starts to collapse and eventually closes in a 
single point from which a thin, needle-like jet is ejected. Existing numerical treatments 
of liquid impact either consider the surrounding air as an incompressible fluid or neglect 
air effects altogether. In contrast, our approach couples a boundary-integral method 
for the liquid with a Roe scheme for the gas domain and is thus able to handle the 
fully compressible gas stream that is pushed out of the collapsing impact cavity. Taking 
into account air compressibility is crucial, since, as we show in this work, the impact 
crater collapses so violently that the air flow through the cavity neck attains supersonic 
velocities already at cavity diameters larger than 1 mm. Our computational results are 
validated through corresponding experimental data. 

Keywords Boundary-integral methods, Roe scheme, solid-liquid impact, supersonic 
air flow 



1 Introduction 

The thin jet ejected after the impact of an object on a liquid surface has been one 
of the icons of fluid mechanics since the days of Worthington more than a century 
ago pQ. Since then a fair amount of computational studies on the impact of liq- 
uid drops [2}{2] or solid objects [5til8j has been reported. In these works a number 
of different methods have been employed including Arbitrary-Langrangian-Eulerian, 
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Fig. 1 The sequence of events as a circular disc of 2 cm radius impacts a water surface at 
1 m/s: (a) Immediately after impact a crown splash is thrown up into the air and an impact 
cavity forms below the free surface. The moving disc draws air into the cavity, (b) Hydrostatic 
pressure pushes the cavity together and the direction of air flow reverses, (c) Eventually the 
cavity closes in a single point, (d) After closure two violent jets are ejected up- and downwards 
from the closure location. The blue and red lines represent the free surface and the disc, 
respectively, from the present numerical simulation with air in (a)-(c). For the jetting in (d) 
the single-phase simulation from 15 is shown. 

Boundary-Integral, Volume-of-Fluid, Level-Set, or combinations thereof. Despite this 
variety a common feature is that the dynamics of the surrounding air was either ne- 
glected altogether 2,3,8-16 or it was treated as an incompressible fluid 4-7,17 in 
the regime for low Mach numbers. 

In this work we report a seemingly simple and harmless situation which nevertheless 
requires to model the dynamics of the gas phase in a fully compressible way: a circular 
disc of 2 cm radius which impacts on a liquid surface with a speed of 1 m/s. Figure [T] 
illustrates the sequence of events during the disc impact extracted from high-speed 
video images 16,18 and compared to the results of our simulations. Upon impact, 
first a thin liquid splash is thrown up all around the circumference of the penetrating 
disc. In the wake of the impactor a large cavity is created which subsequently starts 
to collapse due to the hydrostatic pressure of the surrounding liquid. When the cavity 
closes about half-way down its length two very fast and thin jets are observed shooting 
up- and down from the closure point |15U19tl2T| . 

In the beginning of the process, obviously, air is drawn into the cavity by the moving 
disc. At a later stage, however, this inward flow is counteracted by the shrinking of 
the cavity volume itself and the direction of air flow is not a priori clear. We will show 
that in the competition between cavity expansion (just above the disc) and shrinking 
(around the neck), eventually the shrinking becomes dominant. Accordingly, the air 
flow through the neck reverses and air is pushed out of the cavity. The collapse is so 
violent that the air stream can attain supersonic speeds. To handle this situation, we 
implement an axisymmetric boundary-integral method (BIM) to simulate the motion 
of the liquid surface which is two-way coupled with a fully compressible Roe solver for 
the highly unsteady gas flow. 

In Section [5] we will introduce briefly our BIM and Roe implementations and then 
describe in some detail the two-way coupling between the BIM for the liquid and the 
Roe scheme for the gas domain. Sections 13.11 and 13.21 show the main physical results 
extracted from our simulations and Section 13.31 briefly summarizes the experimental 
validation carried out in [18] . Section [4] concludes the article. 
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2 Numerical methods 

In the impact process viscous effects are negligible as can be seen by estimating the 
Reynolds numbers for gas and liquid as Re g ; = VqRq/ v g j. Here, Vq — 1 m/s is 
the impact speed, Rq = 2 cm is the disc radius, and v g = 1.46 ■ 10~ 5 m 2 /s and 
v\ = 1.12- 10~ 6 m 2 /s are the dynamic viscosities of gas and liquid at 15°C, respectively. 
Both Re s and Re; are larger than 10 3 and viscosity is thus negligible. Furthermore, we 
showed in earlier works [TS] that only a negligible amount of vorticity is present in the 
system. We can thus assume the flow to be inviscid and irrotational as is required for 
the boundary-integral method and the applicability of the inviscid Euler equations. 

Our simulation is split in two stages. In the first, incompressible stage, a two- 
fluid boundary integral method is used to simulate the first part of the impact and 
cavity collapse where gas velocities are of the order of the disc velocity. In the second, 
compressible stage, we use a single-fluid BIM for the liquid coupled to a Roe scheme 
to solve the compressible Euler equations in the gas domain [22] . 

The transition between both stages is fixed at the moment that the air flow through 
the cavity neck reverses and the gas axial velocity at the neck is zero, i.e. n n cck — 0. El 

At the transition moment the liquid cavity has an elongated shape with a neck 
located roughly at the middle as can be seen in Fig. [T](b). Furthermore, the gas flow 
is to a very good approximation one-dimensional directed along the symmetry axis of 
the cavity (as can be verified by flow profiles obtained from the two-fluid BIM shown 
in Section [3TTJ . This makes the situation reminiscent of gas flow through a converging- 
diverging de-Laval nozzle frequently encountered in aerodynamics [23.. In this spirit, 
we remove the inner potential fluid during the compressible stage and replace it with a 
compressible gas described by the one-dimensional Euler equations. To integrate Euler's 
equations in time we use the well-known scheme due to Roe 24 . This description is 
valid along most of the cavity. Above the initial free surface as well as in a small zone 
above the disc the air dynamics can be neglected as will be described in Section T2.3I 
The pronounced difference to the standard nozzle situation, however, is that in our case 
the "nozzle" geometry is determined by a liquid interface which is changing rapidly in 
time and thus creates a highly unsteady gas flow. 

The two-way coupling between the gas and the liquid domains is accomplished via 
(i) the interfacial shape and its instantaneous velocity which is provided by the BIM 
and serves as an input into the gas solver and (ii) the gas pressure which is obtained 
from the solution of the Euler equations and serves as a boundary condition for the 
BIM. Above the location of the initial water level the surface pressure of the BIM 
remains atmospheric. Our combined BIM/Euler method has the advantage that we 
retain the computational efficiency of the BIM which requires only the modeling of the 
boundary and not of the entire liquid domain in contrast to, e.g. compressible air/water 
treatments based on the Volume-of-Fluid method for sloshing tanks [251126] . 

All quantities are non-dimensionalized with the disc radius Kq — 2 cm and the 
impact velocity Vo = 1 m/s. As a third quantity for non-dimensionalization we use the 
density of water p; = 998.23 kg/m 3 for the equations concerning the liquid domain 
and the density of air (at rest under atmospheric pressure) p g = 1.2 kg/m 3 for the gas 
equations. Since viscosity is neglected we have four dimensionless parameters which 



1 The transition point can of course also be taken somewhat later as long as gas velocities 
are still low enough to neglect compressibility. We tried w ncc k — 10, 20, and 50 m/s which all 
give similar results. 
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are the Froude number, the Weber number, the Euler number for the liquid, and the 
Euler number for the gas defined as: 

V 2 

Fr = -IS- (f ) 

gRo 1 ; 

We = (2) 
a 

with the atmospheric pressure p a = 1 Of .3 kPa, the acceleration of gravity g = 9.81 m/s 2 , 
and the air/water surface tension a = 72.8 mN/m. We use cylindrical coordinates r 
and z with z — at the height of the initial free surface as is appropriate for our 
axisymmetric setup. In all computations the disc is assumed to have zero thickness. 

The boundary- integral method for a single fluid (liquid) is briefly sketched in l2.1.1l 
and its extension to two fluids (liquid and air) is given in l2.1.2l Some important aspects 
of our specific BIM implementation are described in 12.1.31 The Roe solver is, again 
briefly, presented in Section [2.21 In Section [2.31 we describe the coupling between the 
gas and liquid domains during the compressible stage. 



2.1 Boundary integral formulation 

2.1.1 Boundary integral method for a single fluid 

If liquid flow is inviscid and irrotational as is the case in our setup 15 , 16 the flow field 
v can be described as the gradient of a scalar potential (j> 

v = V</> (5) 
which satisfies Laplace's equation throughout the liquid domain 

A(j> = 0. (6) 
From Eq. ([6]) the boundary-integral equation can be derived using Green's identities 



,(r) = 



dS' (7) 



with S denoting the boundary, n the normal vector pointing out of the liquid domain, 
and 4> n = n- V</!> the normal derivative of the potential along the boundary. Equation (J7J 
expresses the potential <j> at an arbitrary point r inside the domain (then ft = 4-7r) or 
on its edge (then /3 = 2ir) merely in terms of quantities that are defined on the surface 
of the liquid domain. This has the major advantage that it is sufficient to evolve the 
surface quantities <j> and <j> n effectively reducing the computational problem by one 
spatial dimension. 
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If on every point along the interface either <f> or (f> n is known, Eq. (0 can be solved 
for the missing quantity such that afterwards both (j> and 4>n are available along the 
entire boundary. With this the interfacial velocity can be computed as 

v = + 0„n (8) 

with d/ds denoting the tangential derivative and t the tangential vector along the 
surface. The numerical procedure for this is sufficiently well described in the literature 
[27H29] and the details of the present implementation will be given in Section [2. 1.31 

The boundary conditions for solving Eq. (JT)) are provided as follows: on the disc 
the liquid is required to follow the disc's motion meaning that <j> n = —1 in non- 
dimensional coordinates. On the air/liquid interface the potential <j> is specified by 
integrating Bernoulli's equation: 

rlH'-^M-i-i, (9 ) 

with D/Dt = d/dt+v- V denoting the material derivative, C = Vn the local curvature 
of the interface, and p = Pg/pa the dimensionless gas pressure. Here we assume that 
the pressure at a point far away from the symmetry axis where the flow is quiescent is 
atmospheric, i.e. 1 in non-dimensional coordinates. 

In total, the boundary-integral simulation for a single fluid contains three substeps 
to advance from time step j to j + 1 [271429] . First, with the velocity v^' we integrate 
Bernoulli's equation to compute the potential <^ +1 ^ along the free surface. The 
second step updates the position of the free surface by integrating 

and moves the disc downwards with its prescribed velocity. Finally, we solve the 
boundary-integral equation Q to obtain the liquid potential over the disc <^jj ^ and 

the normal derivative of the potential over the free surface <f)^~^ 1 \ Then the process 
repeats. 



2.1.2 Boundary-integral method for two fluids 

The BIM can be extended to describe two immiscible fluids with a moving interface 
separating both phases. Our approach closely follows that of .'SO TT] and is thus only 
briefly sketched here. The liquid and gas phase both satisfy the boundary-integral 
equation Q within their respective domains. At the interface between the two fluids 
tangential stresses vanish since the fluids are inviscid, i.e., the fluids can slip freely along 
the boundary. To ensure continuity of the interface, however, the normal velocities of 
both phases must exactly balance 

<j>n = ~<t>n,g (11) 

with 4> n ,g being the normal derivative of the gas potential and (f> n that of the liquid as 
defined above. The minus sign is due to the normal vector pointing always out of the 
respective domains. Furthermore, the pressure jump across the interface is given by 
the Laplace pressure. These two conditions are sufficient to derive the two-fluid version 
of the BIM [301I3T] . 
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2.1.3 Specific implementation of the boundary-integral method 

In our axisymmetric situation the surface integrals in Eq. ([7]) can be reduced to one- 
dimensional line integrals in the (r, z)-plane after analytical integration over the az- 
imuthal angle. Numerical integration is carried out using 8-point Gaussian quadrature 
with the weak logarithmic singularities removed analytically as in [29]. Between the 
nodes the interface shape and potentials are interpolated using cubic splines with the 
node-to-node distance serving as the spline parameter s. Natural boundary conditions 
(i.e. a vanishing second derivative with respect to s) are used for all splines, except 
for the potential at the connection between the disc's edge and the free surface as de- 
scribed below. Once the known integrals in Eq. (0 are evaluated we transform Eq. (J7J 
into a matrix equation which is solved by LU decomposition 28,29 . 

Time-stepping for the integration of Eqs. (J9]l and (|10|) is carried out by an iterative 
Crank-Nicholson procedure during the incompressible stage. In the compressible stage, 
we use a simpler forward-Euler scheme for ease of coupling between the BIM and the 
Roe solver in the two respective domains. The time step in the incompressible stage 
is determined by the condition that neighboring nodes may not collide even if their 
velocities were directed exactly towards each other. This leads to: 

At' = f ■ mini (di/vi) (12) 

with i running over all nodes, Vi the free surface velocity at node i, and d.j the distance 
to the neighboring node. This quantity is multiplied with a safety factor / which is in 
most simulations chosen to be 0.1. In the compressible stage the time-step is determined 
by the stability condition of the Roe solver as described in the next section. 

To ensure a continuous boundary of the liquid domain the last node of the free 
surface remains fixed at the disc's edge. This connection point requires some special 
consideration. First, the surface is not smooth and thus the prefactor /3 in Eq. ([7| must 
be modified to become /3 = 2a for the liquid and /3 = 2n — 2a for the gas equations. 
Here, a is the angle connecting the horizontal disc and the tangent to the free surface 
at the connection point through the liquid domain (i.e. a > n in our situation). Second, 
the connection point is at the same time the last node, T, of the free interface and 
the first node, T>, of the disc. To solve Eq. (J7)) we consider it part of the disc, i.e. we 
impose <f> n — — 1 and obtain the corresponding value for (f>. Together with the pinning 
at the disc's edge, the position r and the liquid velocity v at P are thus completely 
determined. We then need to ensure that J- and T> remain identical. For this we first 
copy the spatial coordinates of T> on T . To ensure further that also the velocity v is 
identical on both nodes we project v on the tangential and normal vectors of the free 
surface at T which determines the values of cf> n and d(j>/ds. The latter is imposed as 
a boundary condition on the spline function for (j>. This procedure ensures that the 
velocity v of the connection point is identical when seen from the disc or from the free 
surface, even though the respective normal and tangent vectors are discontinuous. 

In our simulations the free surface extends from the edge of the disc at r = 1 out 
to 100 where any motion is negligible and the surface is cut off. Note that BIMs do not 
require a closed liquid domain since the portion of S at infinity gives no contribution to 
Eq. (|7|) provided that cj> goes to zero there. We use an adaptive mesh to ensure that the 
sensitive areas such as the crown splash or the cavity neck are properly resolved without 
wasting computation time by placing a large amount of nodes on unimportant parts. 
During the incompressible stage the local node distance d is inversely proportional 
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Fig. 2 (a) Illustration of the regridding procedure: every n time steps the surface is recon- 
structed with the new nodes (red crosses) shifted to lie halfway between the old nodes (blue 
circles) in order to avoid amplification of small numerical disturbances, (b) The thin liquid 
sheet ejected after impact breaks up into droplets which are cut off and discarded as illustrated 
by the blue line (circles) before and the red line (crosses) after cut-off. 



to the local curvature C with a proportionality constant between 0.05 and 0.005. We 
impose a maximum distance (dmax = 10) and minimum distance (d m ; n = 0.01) and 
construct our regridding algorithm such that large gradients in the node density, which 
might cause instabilities, are avoided. When coupling with the Euler solver during the 
compressible stage, the BI mesh corresponds to the grid cells used for the Roe scheme 
as will be described in Section 12.31 Note that then we also allow for node distances 
smaller than d m ; n . 

Boundary-integral methods are known to be vulnerable to instabilities (see e.g. [2]) 
due to the lack of a naturally damping viscosity which prevents small numerical dis- 
turbances from building-up over time. To handle such instabilities we use a smoothing 
algorithm as follows: At every nth (usually 2 < n < 10) time step the free surface 
nodes are redistributed such that new nodes fall exactly half-way between old nodes 
as illustrated in Fig. [5] (a). This periodic regridding procedure efficiently ensures the 
stability of the numerical scheme [2J. 

A specific detail of our physical problem is the thin sheet of liquid thrown up around 
the rim of the impacting disc (see Fig. [T](a)). In reality, this sheet will quickly develop 
non-axisymmetric instabilities leading to the formation of individual droplets earning 
it the title of a "crown splash" as seen, e.g., in the famous pictures of [T|. As here we 
are not interested in the details of this splash and our axisymmetric code is not able 
to handle the droplet formation in any case, we cut off the splash as soon the distance 
between the two sides of the liquid sheet at a given point falls below the local node 
distance as illustrated in Fig. [2J (b) . The actual surface surgery is similar to the one 
used in [15] to handle the pinch-off of the cavity prior to jet formation. 

2.2 Roe method for a compressible gas 

In the compressible stage the inner gas is described by the one-dimensional Euler 
equations for conservation of mass, momentum, and energy: 



d (puS) 

at 



+ Eu s 



d{ P S) | d(pus) 

dt dz 

8{pS) d (pu 2 S) 
dz dz 








(13) 



(14) 
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d(pES) d{puHS) „ dS 

All quantities are dimensionless and S is the cross-sectional area of the cavity, p the 
gas density, u the velocity, and p the gas pressure. The total energy E per unit mass 
is defined as 

and the total enthalpy H, again per unit mass, is 

# = (17) 
P 

The fact that we use p a as a reference pressure for p leads to the appearance of Eu g 
in Eqs. (|13[1 - (|15[) . Note the two source terms pdS/dz and —pdS/dt on the right-hand 
side which account for a cavity radius which is changing in space and time. 

Integration of Eqs. (|13p - (|15p is carried out using a Roe scheme [231124] whose 
implementation is fairly standard and thus omitted here. 

The time step during the compressible stage is restricted by the Courant-Friedrich- 
Lewy (CFL) stability condition for the Roe solver. We fix a constant CFL number 
C* < 1 (usually C* — 0.5) and determine the time step by: 

At = C* r (18) 

maxi (|Uj| + Cj) 

with the cell size Az, the local speed of sound in each cell Cj, and the index i running 
over all cells. 

One detail that is of interest are the boundary conditions at the upper and lower 
end of the cavity. In contrast to the standard problem of flow through a nozzle which 
possesses an inlet on one side and an outlet on the other side, in our case gas leaves the 
shrinking cavity on both sides. Since both outflows are subsonic we can prescribe one 
flow quantity at each boundary [23]. To compute the flux through the lower (upper) 
face of the first (last) computational cell we add a boundary cell at each end whose 
state values are calculated as follows. 

At the upper exit we impose that the pressure in the boundary cell be atmospheric. 
In order to compute the density and velocity of the boundary cell, let p^, un, and pj\r 
be the components of the state vector, and cjy = \/JPn/pn the local speed of sound 
in the last computational cell. Assuming that the discharge proceeds isentropically and 
that the 7 + characteristic 

/+ = —?— c + u (19) 

7-1 

transporting information out of the computational domain is conserved, we have for 
the values of the upper boundary (ub) cell: 

l/ 7 

Pub 



Pub^PNi^ (20) 

\pnJ 

2 

u ub = 7 ( C JV - Cub) + U N (21) 

7-1 

Pub = 1- (22) 



At the bottom exit we impose the velocity u e with which the gas leaves the 
computational (Euler) domain. Since velocities between this point and the disc are 
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small (see Section 12. 3[) conservation of mass allows us to calculate this velocity from 
the rate of change in cavity volume between the disc and the boundary cell, i.e. 
u e ■ tv ■ r-g = f„ (pndSh with the surface Si, including the disc itself. Similarly as at 

J Of, 

the upper exit, the state of the lower boundary (lb) cell can then be computed from 
the state of the first computational cell (pi, ui, pi) using conservation of the I~ char- 
acteristic: 



Plb = 7-j" 

c Ib 



Ulb = U e 
Plb = 



7 \ Pi 
Pi 



1/(1-7) 



C/6 = Cl - (Ml - u e ) ■ 



(23) 
(24) 

(25) 
(26) 



These boundary conditions cannot completely prevent reflection of waves at the upper 
and lower ends which causes some small oscillations in the state variables p, u, and p 
as can be seen for example in Fig. [4] (b). The oscillations however remain sufficiently 
small so that the averaged evolution of the gas dynamics can still reliably be extracted. 



2.3 Coupling between boundary-integral and Roe method 

During the compressible stage of the simulation, the gas domain is split into three sep- 
arate zones as illustrated in Fig. [3] (a). The gas pressure which is required for coupling 
to the BIM is obtained in a different way for each zone. First, in the "atmospheric 
zone" the pressure is taken to be atmospheric since the gas dynamics is negligible. 
Next, in the "Euler zone" the pressure is provided by the solution of the Euler equa- 
tions (|13[1 - (|15[) . Finally, in the "bubble zone" the gas dynamics is again neglected. The 
pressure, however, cannot be taken to be atmospheric since the zone is not open to the 
atmosphere. Due to the low gas velocities in that zone (of order of the disc speed) the 
pressure is taken equal to the pressure at the bottom end of the Euler zone. 

The upper end of the Euler zone can remain at a fixed vertical position (between 
and 1 disc radii below the initial free surface). The bottom end is fixed to be at the 
maximum radial extension of the cavity below the neck as depicted in Fig. [3] (a) . Note 
that the Euler zone cannot be extended all the way down to the disc since the free 
surface departs almost horizontally from the disc's edge. This leads to large gradients 
in the cavity radius on the right-hand side of Eq. (| 14p and can thus cause numerical 
instabilities. Since the bottom end of the Euler zone is moving downward in time, its 
length needs to be extended during the simulation. Our algorithm adds a new Roe cell 
whenever the distance between the desired location (maximum radial extension of the 
cavity) and the actual position of the last Roe cell becomes larger than the size of a 
single cell. 

For the results presented here, the initial number of cells in the Roe solver is 600 
and grows dynamically by extension of the Euler zone. To ensure sufficient resolution 
for the high air speeds even during the last moments of the simulation, we double the 
number of cells by splitting each cell in half whenever the minimum cavity radius 
falls below a certain value. This is done twice: at r ncc ^ = 0.2 and r ncc k = 0.05. Due to 
extension and node doubling the number of Roe cells at the end of the simulation is 
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Fig. 3 (a) Illustration of the three zones in which the gas domain is split during the com- 
pressible stage: in the "atmospheric zone" on top, the pressure is always atmospheric and the 
gas dynamics are neglected. In the "Euler zone" we use the Roe scheme to calculate the fully 
compressible gas dynamics. In the "bubble zone" again the pressure is constant and given by 
the pressure at the end of the Euler zone, (b) Schematic illustration of the alignment of Euler 
cells and BI nodes: the BI nodes are always placed at the height of the center of the Euler 
cells. In the neck area a BI node is placed in every 2nd Roe cell, while away from the neck the 
spacing is larger to save computation time. 

somewhat above 3000. The height of the Roe cells is always constant with Az = 0.0012. 
We find that the number of Roe cells is not crucial for the total computation time which 
is mainly determined by the BIM calculations. 

It is crucial to properly align the positions of the BI nodes with the Euler cells. For 
this we place a BI node always exactly in the center of every nth Euler cell as illustrated 
in Fig. [3] (b). Usually n = 2 in a refined zone around the cavity neck and n=5 outside 
this zone. The cross-sectional area S as well as its spatial and temporal derivatives for 
each cell are calculated at the height of the cell center from the splines interpolating 
the cavity surface in the BIM. To avoid numerical instabilities of the BIM we use the 
periodic regridding described above which now makes the BI nodes "jump" between 
Euler cells. This implies that there must always be at least one Roe cell without a BI 
node between two cells which contain a node, i.e. n > 2. 

The two-way coupling between the gas and the liquid domain is accomplished as 
follows. At each time step j we first do a BI step to advance the shape and velocity 
potential of the free surface from j to j + 1 using the gas pressure of step j. The 
BI step further provides the derivatives dS/dz — 2nrdr/dz and dS/dt = 2nrr. This 
calculation is followed by a Roe step using the new cavity shape j + 1 to obtain the 
new pressure at j + 1 and so forth. 

We performed an extensive set of simulations to verify that the results presented 
in the next section are numerically robust when changing any of the above mentioned 
simulation parameters. 

3 Results 

3.1 Justification of the ID compressible scheme 

In Fig. [4] (a) we show the gas velocity through the cavity neck as a function of the 
shrinking cavity neck for the incompressible two- fluid BIM. We use r nec i; instead 

of time to allow for an easier comparison with experiments in Section 13.31 and [IS] . 
Already at r- ncc k = 0.05 (corresponding to 1 mm) an incompressible gas would surpass 
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Fig. 4 (a) The gas velocity at the neck as obtained from an incompressible two-fluid BI sim- 
ulation. The velocity easily surpasses the speed of sound demonstrating the need for our more 
sophisticated compressible approach, (b) The velocity obtained from the multiscale simulation 
(red line) agrees so well with the velocity from the two-fluid BIM (blue line) that both are 
hardly distinguishable for low velocities where compressibility is negligible. The start of the 
red curve marks the transition from the incompressible to the compressible stage. The small 
oscillations in the red curve are due to wave reflections at the ends of the Euler zone, (c) The 
gas flow field obtained from the two-fluid BIM at the moment when the velocity at the neck 
reverses and the simulation passes from the incompressible to compressible stage. Except for 
a rather narrow zone around the neck the flow is to a good approximation one-dimensional. 



the speed of sound. This clearly demonstrates the need for a CFD method which takes 
the fully compressible gas dynamics into account if one wants to study the ejected air 
stream close to cavity collapse. 

A full two-way coupling is required as both gas and liquid flows occur on similar 
time scales: We first estimate the typical time scale for the gas flow T gas as the length 
of the cavity L ~ 10 cm divided by the speed with which a perturbation travels, 
c = 330 m/s, to obtain T gas = 0.3 ms. A typical time scale for the variation of the 
cavity radius can be derived by considering the time that it takes the cavity to collapse 
from a neck radius of 4 mm@ down to zero which is T ca v ~ 1 ms and thus of the same 
order as 1 gas ■ 

For some time after switching from the incompressible to the compressible stage the 
gas velocities are still moderate and compressibility effects should be negligible. We can 
thus expect that in the beginning of the compressible stage the gas velocity obtained 
from our multiscale approach should be similar to the two-fluid BIM in Fig. [4] (a). That 
this is indeed the case is demonstrated in Fig. [5] (b) which gives us a first indication 
that the coupling between the Roe solver and the BIM works correctly. It further gives 
good evidence that the assumption of one-dimensional gas flow in the compressible 
stage is justified. 

To verify the ID assumption more explicitly, Fig.|4](c) shows the flow field obtained 
from the two-fluid BIM at the moment of flow reversal. Except for a small region around 
the neck where by definition the vertical velocity is zero, our assumption is well justified. 



2 At Tneck ?3 4 mm the vertical neck motion starts to reverse which marks the beginning 
when air effects become important, see Fig. 5 of |18| . 
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3.2 Structure of the compressible gas flow 

The intricate structure of the gas flow in the Euler zone is illustrated by the velocity 
profile in Fig. [5] for various instants of time. Here we normalize velocities with the 
speed of sound to obtain the Mach number Ma=u/c. At early times as in Fig. [5] (a) 
one appreciates that at the bottom end of the Euler zone air is entrained by the 
downward moving disc at velocities of the order of the disc speed. This downflux is 
however overcompensated by the shrinking of the cavity around the neck so that the 
total flux is directed upwards as can be seen by the velocity maximum at z ~ —1.9. 
Above the maximum the cavity widens again and the velocity decays towards the upper 
end of the cavity. The consequence of the competition between cavity expansion at the 
bottom and cavity collapse around the neck is the creation of a stagnation point with 
Ma = as is clearly visible in Fig. [5] (a). 

As the collapse progresses gas is pushed through the rapidly diminishing cavity 
neck at ever higher and higher speeds leading to a sharp velocity peak at the neck 
as illustrated in Fig. [S] (b). At the top and bottom boundaries of the Euler zone the 
velocities remain almost unaltered as compared to Fig. 0(a). 

Finally, as the flow speed increases even further a shock wave develops upstream 
of the neck as shown by the magnification in Fig. [5] Thanks to the shock-capturing 
ability of the employed Roe scheme our method is able to handle the shock formation 
quite well. 

Figure [7] shows the Mach number at the cavity neck. For our standard configuration 
of a 2 cm disc impacting at 1 m/s the flow becomes sonic at a neck radius of 0.025 
(corresponding to 0.5 mm). Here we also add data for a higher impact speed of 2 m/s 
which qualitatively shows the same behavior but where sonic flow is attained already 
at a neck radius of 0.06 (1.2 mm). Once sonic velocities at the neck are reached our 
numerical scheme becomes unstable which is why do not present any data beyond this 
point. Note that due to our unsteady situation supersonic speeds are attained even 
earlier at locations above the neck, cf. Fig. [U 

We now turn to study the pressure distribution in the Euler zone which is illustrated 
in Fig. [8] In the early stages of the process (Fig. [8] (a)) the pressure remains virtually 
atmospheric with only a very slight dip around the cavity neck. At a later time, however, 
the pressure at the neck diminishes substantially as can be seen in Fig.[8](b). In a steady 
state situation one would expect the neck pressure to reach a minimum value of 



with 7 = 1.4 the isentropic exponent, as Ma nec k becomes unity. As shown in Fig. (a) 
our situation - although highly unsteady - exhibits a similar behavior with p ncc k ~ 0.6 
at the final instant before the simulation destabilizes. 

Notably, below the neck the pressure is very uniform all the way down to the end 
of the Euler zone (Fig. (a)) which allows us to define a single pressure value for 
the "bubble" between the neck and the disc. Figure [9] (b) demonstrates that in our 
unique situation of air being pushed through a rapidly shrinking nozzle sonic speeds 
can be attained with a bubble pressure which is merely 2% higher than the surrounding 
atmosphere. 




0.53 



(27) 
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Fig. 5 The velocity of the gas stream in the Euler zone for various instants of the cavity 
collapse (right images) with corresponding cavity profiles (left images), (a) In the early stage 
for neck radii around half a disc radius (r^ck = 0.46) the velocity peak is still rather broad. 
At the bottom end gas leaves the Euler zone with a velocity approximately equal to the disc 
velocity of -1 which corresponds to Ma=-0.003. Note that the velocity peak is located somewhat 
upstream of the neck which is marked by the blue dashed line. This is markedly different from 
steady subsonic flow through a fixed nozzle where both would coincide. The location of the 
stagnation point is indicated by the red dash-dotted line, (b) At a later time (r nec i< = 0.05) 
the velocity peak sharpens and increases in magnitude. The peak is now located almost at the 
neck. 




0.5 1 1.5 2 
Ma 



Fig. 6 A shock wave develops (the numbers 1-5 correspond to neck radii of 0.044, 0.036, 
0.034, 0.032 and 0.027; the blue dashed line indicates the neck position for curve number 5). 
For even smaller neck radii the simulation destabilizes and no reliable data is available. 
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Fig. 7 The Mach number at the neck for two different impact speeds: 1 m/s (red solid curve) 
and 2 m/s (green dashed curve). The flow becomes sonic when neck has shrunk to 0.025 or 
0.06 disc radii, respectively. 
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Fig. 8 The pressure profile in the Euler zone at the same instants as in Fig. \E\ The low 
pressure at the neck is caused by the high gas speeds in that region. 



3.3 Experimental validation 

Here we describe briefly the validation of our numerical results with the experimental 
data of [18] which is achieved in three different ways: 

(i) We use smoke particles to directly measure the air speed as it is pushed out of 
the collapsing cavity. We can reliably measure air speeds up to 10 m/s and very good 
agreement with the numerical data is found as presented in Fig. 2 of [18] . 

(ii) In order to confirm the validity of our predictions also for higher gas velocities 
we compare the numerical cavity shape close to pinch-off with experimental images. 
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Fig. 9 (a) The pressure at the neck diminishes over time to reach a minimum value of ap- 
prox. 0.6 atmospheres, (b) In contrast, the pressure deep inside the bubble rises only slightly 
above atmospheric towards the end. 



We find that the experimental shape is not smoothly curved but exhibits a rather 
pronounced "kink" at its neck. This effect - which is due to the low pressure induced 
by the high gas speeds - is not present in single-fluid simulations, but can be reproduced 
rather well by the inclusion of air effects as shown in Fig. 3 of [18] . 

(iii) Finally, we show that the cavity neck in the experiment exhibits a significant 
upward motion prior to final collapse as the fast air stream pushes the surface minimum 
upwards. We find very good agreement for this motion between experiment and our 
multiscale simulations as shown in Fig. 4 of [18| . 

The quantitatively consistent observation of the above air effects in our compress- 
ible simulations and corresponding experiments gives us strong confidence in the reli- 
ability of our numerical scheme. 



4 Conclusions 

We presented a multiscale model to simulate the impact of a solid object onto a liquid 
surface. Our focus was on the fast stream of air that is pushed upwards as the impact 
cavity collapses due to hydrostatic pressure. We showed that in our case of a 2 cm disc 
impacting at 1 m/s the air attains supersonic velocities and thus requires the use of 
a fully compressible computational method - in contrast to existing treatments such 
as Volume-of-Fluid or Level-Set methods [4Tt71ll7] in which the air flow was considered 
incompressible. 

In our simulations the impact process is split in an incompressible and a compress- 
ible stage. During the incompressible stage which covers the first part of the process, 
air is entrained into the cavity at relatively low speeds (compared to the speed of 
sound). This allows us to use a two- fluid boundary-integral method for the gas and the 
liquid domain. The compressible stage starts as the air flow reverses and air is pushed 
out through the cavity neck. In this stage we couple two different methods: a Roe 
scheme to solve the one-dimensional Euler equations in the gas domain and a single- 
fluid boundary-integral method in the liquid domain. The two domains are connected 
via the pressure at the free surface. 

The predominantly one-dimensional character of the air stream and the shape of the 
impact cavity make our system reminiscent to the common problem of air flow through 
a converging-diverging nozzle in aerodynamics. There is, however, an important and 
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fundamental difference: since in our case the confining cavity is a liquid, the "nozzle" 
shape is rapidly evolving in time. 

For low gas velocities we find very good agreement between our simulations and 
direct experimental measurements [TS] . The shape of the cavity as well as a final upward 
motion of the cavity neck [18] give further strong evidence that our multiscale numerical 
method faithfully reflects reality. 
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